# capital-simulacao.R

# Exemplo de uso:
#
# source('capital-simulacao.R')
# source('dados.R')
# SimulaRuina_PremioFixo(500, Danos(5), 250, 4)
# PlotSimulaRuina(500, Anos(5), Danos(5), seq(80,120,length=6), 1.1, 0, 5)
# TempoSimulacao(Danos(5), 100, 4)
#
# PlotCapitalPremio_RuinaFixa(0.1, 0.3, Danos(5), 60, 200)
# PlotPremioCapital_RuinaFixa(.1, .3, Danos(5), 1, 2.5)

source('capital.R')

SimulaRuina_PremioFixo <- function(nsim, X, u, c)
{
	nruinas <- 0
	for(s in 1:nsim) {
		TX <- SimulaGamma(n, X)	
		if(TestaRuina_PremioFixo(TX$T, TX$X, u, c, FALSE))
			nruinas <- nruinas + 1
	}
	nruinas / nsim
}

SimulaRuina_PremioAdj <- function(nsim, X, u, c0, cadj)
{
	nruinas <- 0
	for(s in 1:nsim) {
		TX <- SimulaGamma(n, X)	
		if(TestaRuina_PremioAdj(TX$T, TX$X, u, c0, cadj, FALSE))
			nruinas <- nruinas + 1
	}
	nruinas / nsim
}

ProbSimulaRuina <- function(nsim, X, U, c., cadj)
{
	Y <- rep(length(U), 0)  # result
	TempoSoma <- 0
	TempoNr <- 0
	for(i in 1:length(U)) {
		systime <- unclass(Sys.time())

		if(cadj == 0)
			Y[i] <- SimulaRuina_PremioFixo(nsim, X, U[i], c.)
		else
			Y[i] <- SimulaRuina_PremioAdj(nsim, X, U[i], c., cadj)

		if(i < length(U)) {
			TempoSoma <- TempoSoma + (unclass(Sys.time()) - systime)
			TempoNr <- TempoNr + 1
			print(paste('falta:',
				format((TempoSoma/TempoNr) * (length(U)-i), digits=1), 's'))
		}
	}
	Y
}

## Plot

PlotSimulaRuina <- function(nsim, T, X, U, ccoef, cadj, r)
{
	c. <- Ccoef(T, X, ccoef)
	Y <- ProbSimulaRuina(nsim, X, U, c., cadj)

	plot(U, Y, type='l',
		main=paste('Probabilidade de Ruina, em que',
			'c:', format(c., digits=3),
			', para r:', format(r, digits=5)),
		ylim = c(0, 1),
		cex.main = .9, cex.lab = .8,
		col='blue', lwd=2,
		xlab='Capital', ylab='Probabilidade de Ruina')
}

## Calcular capital/premio necessário para dada probabilidade
## de ruina (usando optimização)

source('progresso.R')

ObtemProbabilidadeRuina_CapitalFixo <- function(prob, Danos, ufixo, cfixo)
{
	oc <- 0
	Steps <- c(200, 10, 1, 0.5, 0.1, 0.05)
	#Steps <- c(10000, 2500)

	for(s in Steps) {
		c. <- oc
		n <- 0
		while(TRUE) {
			c. <- c. + s
print(c.)
			if(ufixo > 0) {
				uarg = ufixo
				carg = c.
			}
			else {
				uarg = c.
				carg = cfixo
			}
			p <- SimulaRuina_PremioFixo(200, Danos, uarg, carg)  # 500
			if(p <= prob)
				break
			oc <- c.
			# se o premio esta a aumentar para valores exorbitantes,
			# passa a frente
			if(s == Steps[1] && c. > 20*Steps[1])
				return(NaN)
		}
	}
	c.
}

PlotCapitalPremio_RuinaFixa <- function(prob1, prob2, Danos, u1, u2)
{
	plen <- 1
	ulen <- 10
	iniciaProgresso(ulen*plen)
	U <- seq(from=u1, to=u2, length=ulen)
	C <- rep(0, ulen)
	media <- mean(Danos)

	P <- seq(prob1, prob2, length=plen)
	Colors <- c(hsv(30/360,1,1), hsv(30/360,1,.8), hsv(30/360,1,.8), hsv(30/360,1,.7), hsv(30/360,1,.6))
	p <- plot
	for(pi in 1:plen) {
		for(i in 1:ulen) {
			C[i] <- ObtemProbabilidadeRuina_CapitalFixo(P[pi], Danos, U[i], 0)/media
			progresso()
		}
		
		p(U, C, type='l',
			main=paste('Ruin probability'),
			cex.main = .9, cex.lab = .8,
			col=Colors[pi], lwd=2,
			xlim=c(u1,u2),
			xlab='Capital (billions in dollars)', ylab='Premium (relative to expected value)')
		p <- lines
	}
	
	lines(c(u1,u2), c(1,1), col='gray')

	legend('topright', format(P, digits=2),
		title='Probability:', title.adj=.05,
		lty=c(1,1), lwd=c(2,2),
		col=Colors, bg='white', cex=.6)
}

## inverso do anterior: premio como objecto; capital como imagem

PlotPremioCapital_RuinaFixa <- function(prob1, prob2, Danos, u1, u2)
{
	plen <- 1
	ulen <- 4
	iniciaProgresso(ulen*plen)
	media <- mean(Danos)
	U <- seq(from=u1, to=u2, length=ulen)
	C <- rep(0, ulen)

	P <- seq(prob1, prob2, length=plen)
	Colors <- c(hsv(30/360,1,1), hsv(30/360,1,.8), hsv(30/360,1,.8), hsv(30/360,1,.7), hsv(30/360,1,.6))
	p <- plot
	
	for(pi in 1:plen) {
		for(i in 1:ulen) {
			#C[i] <- ObtemProbabilidadeRuina_CapitalFixo(P[pi], Danos, 0, U[i])
			progresso()
		}
C[1] <- 23000
C[2] <- 22400
C[3] <- 22100
C[4] <- 20000
print(C)

		if(pi == 1)
			ylim <- c(C[length(C)], C[1])
		p(U, C, type='l',
			main=paste('Ruin probability'),
			cex.main = .9, cex.lab = .8,
			col=Colors[pi], lwd=2,
			xlim=c(u1,u2), ylim=ylim,
			ylab='Capital (in billions of dollars)',
			xlab='Premium (relative to expected value)')
		p <- lines
	}

	lines(c(1,1), ylim, col='gray')

	legend('topright', format(P, digits=2),
		title='Probability:', title.adj=.05,
		lty=c(1,1), lwd=c(2,2),
		col=Colors, bg='white', cex=.6)
}

## Extra: curiosidade

# gráfico com o tempo q demora à medida q se aumenta nro de simulações
# i.e. notação O

TempoSimulacao <- function(X, u, ccoef)
{
	Nsim <- seq(from=10000, to=500, length=6)  # capital
	DT <- rep(length(Nsim), 0)
	c. <- Ccoef(c(1900,2005), X, ccoef)
	for(cadj in c(0.1, 0)) {
		for(i in 1:length(Nsim)) {
			systime <- unclass(Sys.time())
			if(cadj == 0)
				SimulaRuina_PremioFixo(Nsim[i], X, u, c.)
			else
				SimulaRuina_PremioAdj(Nsim[i], X, u, c., cadj)
			DT[i] <- unclass(Sys.time()) - systime
			print(format(DT[i], digits=1))
		}

		if(cadj == 0.1)
			plot(Nsim, DT, type='l',
				main=paste('Estudo da Complexidade das Simulações'),
				cex.main = .9, cex.lab = .8,
				col='blue', lwd=2,
				xlab='Número de Simulações', ylab='Tempo (s)')
		else
			lines(Nsim, DT,
				col='violet', lwd=2)
	}
	
	legend('topright', c('fixed', 'adjusted'),
		title='Premium strategy:', title.adj=.05,
		lty=c(1,1), lwd=c(2,2),
		col=c('violet', 'blue'), bg='white', cex=.6)
}

